#!/bin/sh
#
help(){
echo ''
echo "`basename $0` ver 1.0"
echo ''
echo 'Program to fit the EOS from a VASP run'
echo "Usage: `basename $0` [options] ev.dat [vcol [ecol]"
echo '  vcol and ecol are the column number of volume (Ang^3/atom) and energy (eV/atom), respectively.'
echo '  If not supplied, vcol=1, ecol=2 are assumed.'
echo 'Availables options:'
echo '  -xr xlow xhi  : to define the volume range'
echo ''
exit
}
#
if [ "$#" -lt "1" ]; then
   help
fi
#
noview=0
low='*'
hig='*'
save=0
while [ "$#" -gt "0" ]; do
  if [ "$1" = "-xr" ]; then
     shift; low=$1
     shift; hig=$1
  elif [ "$1" = "-h" ]; then
     help

  elif [ "$1" = "-s" ]; then
     save=1

  elif [ "$1" = "-no" ]; then
     noview=1

  else
     break
  fi
  shift
done
#
vcol=1
ecol=2
if [ $# -ge 2 ]; then
   vcol=$2
fi
if [ $# -ge 3 ]; then
   ecol=$3
fi
#
newf=$1_ev_tmp.dat
cat $1|awk "{if ((\$1+0)==\$1) print \$$vcol, \$$ecol}" > $newf
#
mineng=`sort -n -k 2 $newf|head -1`
E0=`echo $mineng|awk '{print $2}'`
V0=`echo $mineng|awk '{print $1}'`
#
cat > vasp_ev_fit_scpt << EOF
#
EOS(x) = B0*x/Bp * ( (V0/x)**(Bp)/(Bp-1.) + 1.) - B0*V0/(Bp-1.) + E0
#
set fit errorvariables
#
B0 = 1.; Bp = 4.; E0 = $E0; V0 = $V0
fit [$low:$hig] EOS(x) "$newf" u (\$1):(\$2) via B0, Bp, V0, E0

set encoding iso_8859_1
set xlabel 'Average volume (\305^3/atom)'
set ylabel 'Energy (eV/atom)'

set key center top
#set format x '%.1f'
#set format y '%.4f'

set mxtics 5

set term post enha colo 20
set out 'EV_EOS_fit.eps'

plot [$low:$hig] EOS(x) w l lw 2 t 'fitted', "$newf" u (\$1):(\$2) w p pt 7 t 'cal.'

set term png
set out 'EV_EOS_fit.png'
replot

unset out
set term unknown
print ''
print 'B0 = ', B0*160.217733,' +/- ',B0_err*160.217733, ' GPa'
print 'E0 = ', E0,' +/- ', E0_err, ' eV/atom'
print 'B'' = ', Bp, ' +/- ', Bp_err
print 'V0 = ', V0, ' +/- ', V0_err, ' Ang^3/atom'
print ''
print 'Assuming cubic box:'
print 'If  2 atoms per cell, a = ', (V0*2.)**(1./3.), ' Ang'
print 'If  3 atoms per cell, a = ', (V0*3.)**(1./3.), ' Ang'
print 'If  4 atoms per cell, a = ', (V0*4.)**(1./3.), ' Ang'
print 'If  5 atoms per cell, a = ', (V0*5.)**(1./3.), ' Ang'
print 'If  6 atoms per cell, a = ', (V0*6.)**(1./3.), ' Ang'
print 'If  7 atoms per cell, a = ', (V0*7.)**(1./3.), ' Ang'
print 'If  8 atoms per cell, a = ', (V0*8.)**(1./3.), ' Ang'
print 'If 10 atoms per cell, a = ', (V0*10.)**(1./3.), ' Ang'
print 'If 12 atoms per cell, a = ', (V0*12.)**(1./3.), ' Ang'
print 'If 16 atoms per cell, a = ', (V0*16.)**(1./3.), ' Ang'
print ''
EOF
#
gnuplot vasp_ev_fit_scpt
if [ "${noview}" -lt 1 ]; then
   evince EV_EOS_fit.eps
fi

if [ "$save" -eq 0 ]; then
   rm -rf EV_EOS_fit.eps vasp_ev_fit_scpt $newf
fi
rm -rf fit.log
#
exit 0
